In [1]:
%matplotlib inline
#%matplotlib widget
from astropy.cosmology import LambdaCDM
import numpy as np
import matplotlib.pyplot as plt
from astropy import constants as const
import astropy.units as u
from scipy.integrate import quad
import ezgal # BC03 model maker
import os
In [2]:
cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7, Tcmb0=2.725)
In [3]:
# check to make sure we have defined the bpz filter path
if not os.getenv('EZGAL_FILTERS'):
os.environ['EZGAL_FILTERS'] = (f'{os.environ["HOME"]}/Projects/planckClusters/MOSAICpipe/bpz-1.99.3/FILTER/')
model = ezgal.model('bc03_ssp_z_0.02_salp.model')
model = model.make_exponential(1)
model.set_cosmology(Om=cosmo.Om0, Ol=cosmo.Ode0, h=cosmo.h, w=cosmo.w(0))
model.add_filter('g_MOSAICII.res', name='g')
model.add_filter('r_MOSAICII.res', name='r')
model.add_filter('i_MOSAICII.res', name='i')
model.add_filter('z_MOSAICII.res', name='z')
model.add_filter('K_KittPeak.res', name='K')
# Blanton 2003 Normalization
Mr_star = -20.44 + 5 * np.log10(cosmo.h) # abs mag.
# set the normalization
model.set_normalization('sloan_r', 0.1, Mr_star, vega=False)
In [4]:
# desired formation redshift
zf = 6.0
# fetch an array of redshifts out to given formation redshift
zs = model.get_zs(zf)
# Calculate some cosmological stuff
DM = cosmo.distmod(zs)
dlum = cosmo.luminosity_distance(zs)
Need to compute the cluster volume...
$M_{vir} = 4/3 \pi r^3_{vir} \rho_c(r<r_{vir}) = 4/3 \pi r^3_{vir} \Delta_c \rho_c$
if we let $\Delta_c = 200$ then
$M_{200} = 4/3 \pi r^3_{200} 200 \rho_c$ with $\rho_c = \frac{3H(z)^2}{8\pi G}$
or just $M_{200} = V_{200}200\rho_c$. So we'll make a function to calculate $\rho_c$. And we'll make use of the astropy units package to do all the unit analysis for us.
Don't forget that $H(z) = H_0E(z)$
The Schechter Function:
For Luminosity:
$\Phi(L) = \phi^\star \frac{L}{L_\star}^\alpha e^{-\frac{L}{L_\star}}$
For Magnitudes:
$\Phi(M) = \phi^\star\frac{2}{5}log(10) (10^{\frac{2}{5}(M_\star - M)})^{\alpha+1} e^{-10^{\frac{2}{5}(M_\star - M)}}$
In [5]:
def rho_crit(z, cosmo):
# convert G into better units:
G = const.G.to(u.km**2 * u.Mpc/(u.M_sun * u.s**2))
return 3 / (8 * np.pi * G) * cosmo.H0**2 * cosmo.efunc(z)**2 # Mpc^3
def schechterL(luminosity, phiStar, alpha, LStar):
"""Schechter luminosity function."""
LOverLStar = (luminosity/LStar)
return (phiStar/LStar) * LOverLStar**alpha * np.exp(- LOverLStar)
def schechterM(magnitude, phiStar, alpha, MStar):
"""Schechter luminosity function by magnitudes."""
MStarMinM = 0.4 * (MStar - magnitude)
return (0.4 * np.log(10) * phiStar * 10.0**(MStarMinM * (alpha + 1.)) * np.exp(-10.**MStarMinM))
In [6]:
from astropy.table import Table
from scipy.interpolate import interp1d
z1 = 0
z2 = 2
dz = 0.025
# build the mass array
zarr = np.arange(z1, z2 + dz, dz)
ps2 = Table.read('../../catalogs/PSZ2v1.fits')
df2 = ps2.to_pandas()
data = df2[['REDSHIFT', 'MSZ']]
data['REDSHIFT'].replace(-1, np.nan, inplace=True)
# redshift bins
zbins = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 3]
nMasses = 100
big_mass = []
for j in range(nMasses):
mass = np.ones_like(zarr) * 1e14
for i in range(len(zbins) - 1):
mask = (zbins[i] <= zarr) & (zarr < zbins[i + 1])
mass[mask] *= float(data.loc[(zbins[i] <= data['REDSHIFT']) & (data['REDSHIFT'] < zbins[i + 1]), 'MSZ'].sample()) * cosmo.h
big_mass.append(mass)
mass = np.vstack(big_mass)
mass_func = interp1d(zarr, np.median(mass, axis=0))
In [7]:
# So now we are going to calculate the volumes as a function of z
#M200 = mass_func(zarr) * u.solMass
M200 = 1e15 * u.solMass
V200 = M200/ (200 * rho_crit(zs, cosmo))
# Calculate the M_star values
Mstar = model.get_absolute_mags(zf, filters='i', zs=zs)
# calculate the abs mag of our limiting magnitude as a function of z
mlim = 23.5
#Mlim = Mstar - 2.5 * np.log10(0.4)
Mlim = mlim - cosmo.distmod(zs).value - model.get_kcorrects(zf, filters='i', zs=zs)
# Here are the Schechter function stuff from Liu et al.
phi_star = 3.6 * cosmo.efunc(zs)**2
alpha = -1.05 * (1 + zs)**(-2/3)
fr = 0.8*(1 + zs)**(-1/2)
#alpha = np.ones_like(alpha) * -1
#Mpiv = 6e14 * u.solMass
#zpiv = 0.6
#alpha = -0.96 * (M200 / Mpiv)**0.01 * ((1 + zs)/ (1 + zpiv))**-0.94
#phi_star = 1.68 * (M200 / Mpiv)**0.09 * ((1 + zs)/ (1 + zpiv))**0.09 * cosmo.efunc(zs)**2
#fr = 0.62 * (M200 / Mpiv)**0.08 * ((1 + zs)/ (1 + zpiv))** -0.80
In [8]:
LF = []
for phi, a, M_star, M_lim in zip(phi_star, alpha, Mstar, Mlim):
if M_lim < M_star - 2.5 * np.log10(0.4):
Mlimit = M_lim
else:
Mlimit = M_star - 2.5 * np.log10(0.4)
y, err = quad(schechterM, -30, Mlimit, args=(phi, a, M_star))
#print(M_star - M_lim, y)
LF.append(y)
In [9]:
plt.figure()
plt.plot(zs, (LF * V200.value + 1) * fr)
ax = plt.gca()
ax.set_yticks(np.arange(0, 75, 10))
plt.xlim(0.1, 5)
plt.ylim(0, 80)
plt.xlabel('redshift')
plt.ylabel('N (r < r200)')
plt.grid()
In [10]:
# calculate the abs mag of our limiting magnitude as a function of z
mlim = 23.5
#Mlim = model.get_absolute_mags(zf, filters='i', zs=zs) - 2.5 * np.log10(0.4)
Mlim = mlim - cosmo.distmod(zs).value - model.get_kcorrects(zf, filters='i', zs=zs)
plt.figure()
plt.plot(zs, model.get_absolute_mags(zf, filters='i', zs=zs), label='Lstar')
plt.plot(zs, Mlim, label='Mlimit')
plt.plot(zs, model.get_absolute_mags(zf, filters='i', zs=zs) - 2.5 * np.log10(0.4), label='0.4Lstar')
plt.grid()
plt.xlabel('redshift')
plt.ylabel('abs Mag')
plt.legend()
Out[10]:
In [11]:
Mlim
Out[11]:
In [12]:
Mstar - 2.5 * np.log10(0.4) # 0.4L* magnitudes
Out[12]:
In [13]:
np.array(LF) # LF integration output
Out[13]:
In [14]:
alpha
Out[14]:
In [15]:
phi_star
Out[15]:
In [16]:
fr # red fraction
Out[16]:
In [17]:
zs # redshift array
Out[17]:
In [18]:
V200.value # cluster volume
Out[18]:
In [7]:
200 * rho_crit(zs, cosmo)
Out[7]:
In [12]:
plt.plot(zs, (V200/(4/3 * np.pi))**(1/3))
Out[12]:
In [ ]: